Simulation of tunneling in the quantum tomography approach 
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The new method for the simulation of nonstationary quantum processes is proposed. The method 
is based on the tomography representation of quantum mechanics, i.e., the state of the system is 
described by the nonnegative function (quantum tomogram). In the framework of the method one 
uses the ensemble of trajectories in the tomographic space to represent evolution of the system 
(therefore direct calculation of the quantum tomogram is avoided). To illustrate the method we 
consider the problem of nonstationary tunneling of a wave packet. Different characteristics of 
tunneling, such as tunneling time, evolution of spatial and momentum distributions, tunneling 
probability are calculated in the quantum tomography approach. Tunneling of a wave packet of 
composite particle, exciton, is also considered; exciton ionization due to the scattering on the barrier 
is analyzed. 

PACS numbers: 03.65.Wj, 02.70.-Ns, 03.65.Xp 



I. INTRODUCTION 

Nowadays simulation of quantum systems is developed 
to a high extent (see, e.g., reviews 0, 0). However, 
common simulation methods, for example, Path Integral 
Monte Carlo or Wigner dynamics, use non-positively de- 
fined functions (wave function, the Wigner function, etc.) 
to describe a quantum state. This leads to the difhcul- 
ties with convergence of corresponding integrals, espe- 
cially harmful for the simulation of Fermi systems (sign 
problem) . There is a hope that using a real nonnegative 
function, describing the quantum state, one can avoid 
these difficulties. 

The real nonnegative function in the phase space, com- 
pletely describing the quantum state, was proposed 60 
years ago (j^, see also During the last decade 

another very interesting representation has been actively 
developed: the quantum tomography, operating with the 
ensemble of scaled and rotated reference frames, instead 
of the phase space H S i, H El 113 • 

In the framework 

of this formalism the state-describing function (called 
marginal distribution or quantum tomogram) is real and 
nonnegative. The advantage is that the quantum tomo- 
gram is a probability distribution shown to completely de- 
scribe the quantum state 0, . It is one of the reasons 
why the quantum tomography has become so popular. 

In this paper we propose a new method for computer 
simulation of nonstationary quantum processes based on 
the tomography representation of quantum mechanics 
and illustrate it considering the problem of nonstationary 
tunneling of a wave packet. Many simulation approaches 
in nonstationary quantum mechanics are based on the 
numerical solution of the time-dependent Schrodinger 
equation. There are also methods using the ensembles 
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of classical trajectories to simulate quantum evolution. 
For example, the method of "Wigner trajectories", based 
on the Wigner representation is well known (see, 
e.g., review [l| for details) and was recently successfully 
app lied to investigate the tunneling of a wave packet 
Using such trajectories, one does not have to 
store the large arrays representing, say, the wave func- 
tion (contrary to grid methods), therefore there is a gain 
in computer memory. 

The quantum tomogram w depends on the variables 
{X,^,iy}, where X = + vp, and q,p are the coor- 
dinates and momenta of the system, respectively, /z, v 
are the parameters of scaling and rotation of reference 
frame in the phase space. The quantum tomogram is 
nonnegative and normalized in X direction, therefore it 
can be interpreted as a distribution function of the value 
X. In our method the ensemble of trajectories in space 
{X,^,v} is introduced to describe the quantum evolu- 
tion. The trajectories are governed by the dynamical 
equations obtained from the evolution equation for the 
quantum tomogram. 

We demonstrate the method considering the nonsta- 
tionary tunneling of a wave packet through the poten- 
tial barrier. For this problem we calculated tunneling 
times, which are of interest nowadays, both for basic sci- 
ence {e.g., what is the time spent by an atom to tunnel 
from the trap?) and for applications (electronic tunnel- 
ing time is connected with the operation rate of some 
nanostructure-based devices). We also analyzed the evo- 
lution of the wave packet in coordinate and momentum 
spaces in details. Another demonstration was designed 
to show that our method is not restricted to one-particle 
simulations. Namely, we investigated the tunneling of 
a wave packet of composite quasiparticle, an exciton 
(coupled electron and hole in semiconductor), in a one- 
dimensional nanostructure (quantum wire). In this case 
we had two degrees of freedom. For this problem, in 
addition to the probability density evolution, we deter- 
mined the probability of ionization due to electron and 
hole scattering on the barrier in different directions. 
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In Sec. m we present description of the method, in 
Sec, mil describe the model problem and main results for 
a wave packet tunneling. Exciton tunneling is considered 
in Sec. II VI and the work is summarized in Sec. IVI 



II. THE METHOD OF SIMULATION 

The quantum tomogram wiX, u, v ) is connected with 
the density matrix p{q,q') as p^H^: 

piq,q')= j z„(X,^,g-g')e'(^-^(^+«')/2)^,(l) 



e 



i{k{X^ij.q~up)+pu) ( .^^ M dpdkdqdu 



(2) 



Consider the case of the particle with mass m in one- 
dimensional space. If the Hamiltonian of the system is 



H 



2m 



V{q). 



(3) 



then the integral transformation Q applied to the time- 
dependent evolution equation for the density matrix gives 



ii dw ^dV(q) f V d 

.0.. . 1M a„2n+1 1 o— ) ^ = 0' (4) 



-(2n+l)! dq 



2 dX 



where we use h ^ \ and q is given by 



d 
'dX 



d_ 
9/i 



(5) 



Eq.l^J can be rewritten as 



dw dw 
'dt'^dX 



dw 
dw 



+ —^x{X, fi, v) + -7^G^,{X, fi, v) + 
d^i 



di 



G,iX,p,i^)^0, (6) 



where functions G depend on quantum tomogram, its 
derivatives and antiderivatives (the latter corresponding 
to terms with {d/dX)~^ in Eq.Q). Generalization for 
the case of more variables is straightforward because the 
form of the equations does not change. Functions G for 
the problem under investigation are given in Sec. lIIII The 
evolution equation rewritten as lO has the form of con- 
tinuity equation for the quantum tomogram 



dw dw dw ■ 
~ ~dt dX 



dw . dw . ^ 



This equation is analogous to the continuity equation 
for classical distribution function and Liouville equation. 



As known, the characteristics of Liouville equation are 
the classical trajectories in phase space and they obey 
Hamilton equations of motion. The quantum tomogram 
is nonnegative and we use it as a distribution function 
for trajectories in the space {X, /i, v}, obeying the equa- 
tions analogous to Hamilton equations for the classical 
trajectories. From the comparison of Eq.® with Eq.© 
it is obvious that the trajectories are governed by the 
equations 

X = GxiX, J^), A = G^{X, ^i,iy),i^ = G,{X, fi, v) (8) 

The trajectories are used to avoid the direct calcula- 
tions of the distribution function (contrary to grid meth- 
ods where the wave function is calculated at every point 
to solve numerically Schrodinger equation). Hence it is 
necessary to use some approximation for the quantum 
tomogram and we use local exponential approximation 
(as in Ref. 15] for the Wigner function): 

(9) 

where y — {X, /i, ^'}, and ya is the point under con- 
sideration. Parameters of this approximation are ma- 
trix Aa and vector ha, and some combinations of these 
parameters enter the evolution equation instead of 
the derivatives and antiderivatives of the quantum tomo- 
gram. Calculation of average X,fi,v and their average 
products allows to obtain Aa and ba- After that, func- 
tions G are known and dynamical equations © can be 
solved numerically. 

We would like to emphasize that we use the local ap- 
proximation (O only for the calculation of r.h.s. in the 
equations of motion ||SJ|. The use of ensemble of tra- 
jectories to represent the quantum tomogram means the 
approximation of quantum tomogram as a set of delta- 
functions, each delta- function corresponds to one tra- 
jectory. If the number of trajectories approaches infin- 
ity, the quantum tomogram can be approximated by the 
set of delta-functions with arbitrary precision. This is 
analogous to what is conventionally done in classical sta- 
tistical mechanics for the distribution function in phase 
space. But unlike the classical statistical mechanics now 
the trajectories are not independent: the approximation 
© is used to take the non-local character of quantum- 
mechanical evolution into account. 

The validity of this approximation holds if the quan- 
tum tomogram is smooth and the trajectories are close 
to each other. For example, this approximation can fail if 
one tries to consider a plain wave with wave vector k: in 
this case w{X, fj, = 0,v = 1) = S(X — k). Approximation 
^ also works not well for unbounded motion, because 
the trajectories scatter with time. If there arc few tra- 
jectories in the region around the given point, then the 
approximation @ will not reconstruct the quantum to- 
mogram well, due to lack of statistics. 

We consider the tunneling of wave packets through the 
barrier, and comparing our results with exact quantum 
computation we see that approximation (jS)) holds for this 
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problem (see Sees. IIIIIIVII . For example, considering the 
tunneling through the potential barrier from the well, 
we deal with both the region of bounded motion (in 
the well) and unbounded motion (beyond the barrier), 
still the approximation works not bad (Sec. . For 
higher initial energy the penetration through the barrier 
is higher, and the evolution of most trajectories corre- 
sponds to unbounded motion, then the validity of Eq.@ 
becomes poorer (see the end of Sec. in agreement 
with the discussion above. But in general local approx- 
imation works satisfactory even for quite long time 
intervals (see Fig. . 

To obtain any information about the system, we have 
to calculate some average values. Consider an arbitrary 
operator A(q,p). Average value (A) of corresponding 
physical quantity is calculated in the tomographic rep- 
resentation of quantum mechanics as [l9| 

{A) = j A{ji,u)e'^w{X,^ji,v)dXd^idv, (10) 

where A(/x, v) is the Fourier component of the Weyl sym- 
bol A^ {q,p) of operator A{q,p) (see, e.g., |lj): 

A{^l, u)= j A'^{q,p)exp{-i{m + i^p))'^ (11) 

For the calculation of average values we use the follow- 
ing approximation of quantum tomogram: 

J 

w{X, >y,t) = Y^ S{X - X,{t))5{pi - fij{t))6{,^ - ,y,{t)), 

(12) 

where the summation is made over all J trajectories; 
Xj{t), ijij{t),Vj{t) are the coordinates of the j-th trajec- 
tory in {Xjfijv} space at time t. Such approximation 
corresponds to use of the ensemble of trajectories. In the 
regions, where w{X, fi, v) is small, trajectories are rare, 
and where it is large, trajectories are accumulated. The 
more trajectories are used, the better the approximation 
(|12|1 works. If during the simulation the wave function 
has the form of a compact wave packet, even consisting 
of several distinct parts, approximation 1)12(1 holds, be- 
cause in this case one has the compact sets of trajectories 
providing good statistics. This is the case for problems 
considered, and therefore the use of this approximation 
does not change results essentially, in comparison with 
the exact quantum computation (Sees. irTllIVp . 

For the operators A{q), depending on q only, expres- 
sion for {A) takes the form: 

{A) = j A{X)w{X, ^ = 1, = 0)dX, (13) 

where A{X) is the function corresponding to the operator 
A{q) in coordinate representation, A{X) — A{q = X). 
The method of calculation of an average {A{q)) at arbi- 
trary time with the approximation H12|) . is quite sim- 
ple. One just takes into account the trajectories with 



any X and with fj,(t),h'(t) from the small region near 
fi = I, v ^ only, and performs a summation of A{X) 
over all such trajectories. 

The developed method is similar to well-known method 
of Wigner trajectories (see l| for review), where the en- 
semble of trajectories is introduced in the phase space, 
with the Wigner function used as a quasi-distribution 
function. The quantum tomogram is defined in the space 
{X,fi,i>}, that is not as simple for understanding as the 
phase space used in the Wigner approach. On the other 
hand, quantum tomogram is true distribution function, 
while the Wigner function can be both positive and neg- 
ative. In spite of these differences the two approaches 
are quite close in general, and there are some difficulties 
common for both methods. The discussion (see above) of 
the approximation for the tomogram (such as lO) con- 
cerns, in principle, the method of Wigner trajectories 
either. Another important example is the discontinuity 
of the Wigner trajectories discussed in Ref. |23|, that is 
due to the fact that the trajectories are not independent 
as in classical statistical mechanics, their evolution de- 
pends on their local distribution. The same is applicable 
for the quantum tomography approach, where the tra- 
jectories are also not independent for the same reason. 
In Ref. fU a possible alternative to Wigner function was 
proposed. Namely, in that work the authors proposed to 
use the Weyl transforms of some operators (the Wigner 
function, up to the constant, is the Weyl transform of 
the density operator) instead of the Wigner function to 
generate the ensemble of trajectories. The same can be 
introduced for the quantum tomography. Applying the 
transform ^ to the matrix elements of an operator A 
we obtain the symbol wa {X, fi, u) . This function is not 
nonnegative in general, but in analogy with the Wigner 
trajectories one can use wa(,X, ji, v) to develop some new 
ensemble of trajectories. Probably, as with the Weyl 
transforms trajectories 20], it will be more convenient to 
use the trajectories corresponding to wa in certain cases. 
Of course, this problem needs further investigation. 

III. SIMULATION OF TUNNELING OF A 
WAVE PACKET 

A. The model and calculated average values 

We choose the external potential to coincide with the 
potential used in for comparison of the results of 
simulation in quantum tomography approach with those 
obtained by other methods. We consider behavior of the 
wave packet in one-dimensional space in external poten- 
tial 

Vi,) = (14) 

As the potential has only the second and third powers 
of coordinate, all its derivatives of order more than the 
third vanish. Evolution equation in this case has the form 
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(n = l) 



dw fi dw 2 
dt m dv 



dViq) 



1 d^Vjq) fv d 
6 dq3 \2dX 



V d 
2'dX 

3" 



w = 



(15) 



For the potential given by Ea. p4ll evolution equation 



reads as: 



dw /J, dw 



at m ov a/i 



d_ 

dX 



12 dX^ 



(16) 



and dynamical equations have the form 
dX bu^ 1 d^w 



dt 

9/i 

'dt 

dv 

'dt 



12 w dX'^ 
o bv 

w 



d 

dX 



dw 
dfj, 



m 



(17) 



We use atomic units throughout, h = rUe = |e| — 
1, where mg and e are the mass and charge of a free 
electron. The particle with mass m = 2000 is regarded. 
Parameters of the potential are luq — 0.01 and b = 0.2981. 
This potential has the minimum at g = {V{Q) — 0) and 
maximum ai q = 0.6709 (t/(0.6709) = 0.015), therefore 
here we consider the motion of a particle in the potential 
well with infinite left wall and the barrier of height 0.015 
at g = 0.6709. This model problem roughly describes 
nonstationary tunneling of an atom from the trap. 

Initially the particle represented by the wave packet is 
situated to the left from g = 0, its mean momentum is 
zero. The particle can oscillate in the potential well and 
can tunnel or pass above the barrier. The probabilities of 
these processes depend on the initial energy of the wave 
packet. We consider the problem, where all parameters, 
except the initial mean coordinate go of the wave packet, 
are fixed (initial mean momentum equals zero, disper- 
sions of the wave packet in coordinate and momentum 
spaces are « 0.3 and « 1.6, respectively). 

We solve the equations (|17|l numerically. As in ^T^] we 
consider three values of go^ —0.2,-0.3 and —0.4. The 
most interesting quantities characterizing tunneling are 
reaction probability and tunneling time. Reaction prob- 
ability is defined as 



\il'{x, t)\'^dx^ 



(18) 



where qa = 0.6709 (the point where potential has the 
maximum), the maximum value of reaction probability 
is unity. Reaction probability shows what part of the 
wave packet is currently beyond the barrier. 



Tunneling time of the wave packet is also an important 
feature of tunneling. There are a lot of methods to de- 
termine tunneling time HHHEllllllllllllllllIll 

We use the approach where tunnel- 
ing time is calculated as the difference of presence times 
(see [s^ l for review) at point Xa and X},, located on the 
opposite sides of the barrier: 



tT{Xa,Xi,) = {t{xb)) - (t{Xa)) 

The presence time at arbitrary point xq is 



J t\%l}{xo,t)\^dt 



Smxo,t)\Ut 





(19) 



(20) 



B. Reaction probability 

Here we present the results obtained by means of our 
method and compare them with the exact numerical so- 
lution of Schrodinger equation. In Fig. ^ we present 
the dependence of reaction probability H18|) on time for 
three values of initial mean coordinate of the wave packet 
go — —0.2, —0.3 and —0.4, corresponding mean ener- 
gies of the wave packet are w 0.75Vo, ~ 1.25Vo and 
K, 2.0Vb, respectively. Solid lines represent the results of 
simulation in the quantum tomography approach (QT) 
and dashed lines correspond to the numerical solution 
of Schrodinger equation (exact quantum computation). 
Due to the increase of initial mean energy with the in- 
crease of I go I, the portion of high energy components in 
the wave packet grows. This leads to the higher portion 
of components, which pass through the barrier, either be- 
cause their energy is greater than the height of the bar- 
rier, or due to tunneling. Therefore, with the growth of 
go I , reaction probability becomes larger and one can see 
that the curves corresponding to different go are above 
each other in Fig.^ The time evolution of reaction prob- 
ability is qualitatively the same for every go. The com- 
ponents, which have passed through the barrier, can not 
return, because for g > 0.6709 potential diminishes with 
the growth of coordinate, and so reaction probability can 
not decrease with time. At first it grows rapidly due to 
transmission of components with the energy higher than 
the height of the barrier (comparison with the classical 
solution of corresponding problem, for which only trans- 
mission above the barrier is possible, convinces us in it). 
Then reaction probability continues to grow slowly, be- 
cause of the tunneling. All these features present both 
for QT simulation and for exact quantum computation. 

In comparison with the exact computation, reaction 
probability for the QT simulation is slightly higher. Note 
also some difference in the character of increase of the 
reaction probability for QT simulation and exact solu- 
tion: in the former case the curves are not so smooth. 
These differences are due to the finite number of tra- 
jectories used in QT simulation: for smaller number of 
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-1.0 -0.5 0.0 0.5 
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FIG. 1: The dimensionless reaction probabilities l)18|l for three 
values of initial mean coordinate of the wave packet: qo = 
—0.2, —0.3, and —0.4 a.u. Solid lines are for the simulation 
in quantum tomography approach, dashed lines are for the 
exact numerical solution. 



FIG. 2: Probability density in coordinate space for QT sim- 
ulation (histograms) and exact solution (smooth lines), at 
times t = a.u. (left) and t = 200 a.u. (right). The bar- 
rier is at the point 0.6709 a.u., qo = —0.2 a.u. 



trajectories (not shown) reaction probability curves re- 
semble staircase more evidently (this is connected with 
the overestimation of the role of wave packet oscillations 
in the well for the finite number of trajectories), and 
quantitative deviation from exact result is stronger. But 
in general, for quite large number of trajectories, as for 
the case shown in Fig. ^ QT simulation results on reac- 
tion probability are quite close to those obtained through 
the exact quantum computation (compare also with the 
method of "Wigner trajectories" in the work by Donoso 
and Martens p^). 



Evolution of the wave packet and tunneling 
times 
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X (a.u.) 



Besides the reaction probability we also obtained a 
number of new qualitative and quantitative results, which 
described in details the behavior of the wave packet dur- 
ing tunneling. We also calculated tunneling times using 
the concept of presence time (see below). 

The following discussion concerns tunneling of the 
wave packet with initial mean coordinate qo = —0.2. 
Note that the maximum of the potential is at the point 
X = 0.6709. We present the normalized probability den- 
sity 1-0(2;) 1^ in coordinate space (Figs.[5]-0 and |-0(p)p in 
momentum space (Figs.|31,El) fo^' several successive time 
moments. In these figures smooth lines show the shape 
of the wave packet obtained by means of exact quantum 
computation. Histograms represent the result of single 
QT run. One can consider many runs with the same 
number of trajectories and average the probability den- 
sity over all these runs to obtain smoother picture. But 



FIG. 3: Probability density in coordinate space for QT sim- 
ulation (histogram) and exact solution (smooth line), at time 
t — 300 a.u. The barrier is at the point 0.6709 a.u., go = —0.2 
a.u. 



here we would like to show what QT simulation can give 
for one run, in comparison with the exact quantum com- 
putation. Therefore the histograms (QT) fit the smooth 
solid lines in Figs. |21- El (exact solution) not ideally, still 
the resemblance is obvious. 

First, consider the probability density |'!/)(a;)|^ in co- 
ordinate space (Figs. 121-©. One can see (Fig. that 
initially the wave packet has Gaussian form. It begins 
to move as a whole towards the potential minimum at 
X — (initial mean momentum is zero, but the potential 
inclines in that direction), passes that point, is acceler- 
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FIG. 4: Probability density in coordinate space for QT sim- 
ulation (histogram) and exact solution (smooth line), at time 
t — 400 a.u. The barrier is at the point 0.6709 a.u., go = —0.2 
a.u. 



ated, and collides with the barrier. During the motion 
the wave packet broadens (due to dispersion in momen- 
tum space, compare right and left plots in Fig. ^ but 
the interaction with the barrier changes its form more 
substantially {t = 300 and t = 400, Figs.ElEJl. The wave 
packet shrinks a little, some components pass through 
the barrier and transmitted part can be seen beyond the 
barrier (x — 0.6709). As the transmitted part can not 
return and accelerates (potential diminishes with the dis- 
tance for X > 0.6709), the enriching of the wave packet 
by high-energy components must present (see below). 

All features described in the previous paragraph are 
common for both the exact solution and QT simulation. 
The histograms in Figs. |21- El fit the smooth lines repre- 
senting exact solution better for earlier times, but even 
after the interaction with the barrier (Fig. 0)), resem- 
blance is quite close. This shows that approximations 
and H12|) work not bad for the problem under consid- 
eration. 

Proceed now to the evolution of the wave packet in 
momentum space. To confirm our analysis on the ac- 
celeration of the transmitted part of the wave packet, 
we present the probability density \'ip{p)\'^ in momentum 
space in Figs. |S1 at times t = Q and t — 400, respec- 
tively. As the wave function is initially the Gaussian 
wave packet, the initial distributions both in coordinate 
and momentum space are Gaussian (compare Fig. [21 and 
Fig. O . But after the wave packet has interacted with 
the barrier, the distribution in momentum space changes 
substantially (Fig. jSJ. The higher the energy of incident 
particle, the greater the tunneling probability. There- 
fore the barrier transmits mainly the wave packet com- 
ponents with relatively high energy, serving as an energy 
selector. Components with high momentum arise, as the 



0.15 




p (a.u.) 



FIG. 5: Initial probability density in momentum space for 
QT simulation (histogram) and exact solution (smooth line), 
t = a.u. 



0.15 




-10 10 20 30 

p (a.u.) 

FIG. 6: Probability density in momentum space for QT simu- 
lation (histogram) and exact solution (smooth line) at f = 400 
a.u. 



transmitted part of the wave packet is accelerated in the 
region of lowering potential beyond the barrier, and we 
observe the enriching of the wave packet by high energy 
components. The resemblance between the histograms 
(QT simulation) and smooth solid lines (exact solution) 
is somewhat poorer for momentum distribution at large 
times {t — 400, Fig. O than for coordinate distribution 
(Fig. 21 . This is due to the fact that one deals with fi- 
nite number of trajectories and has to sample quite large 
interval in momentum space, because the transmitted 
part is permanently accelerated. Therefore, with time, 
momentum distribution spreads and there are not many 
trajectories with fi, v close enough to = 0, = 1 for 
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qo (a- u) 

FIG. 7: Tunneling times with errors for several values of 
initial mean coordinate of the wave packet qo ■ Results of the 
QT simulation (squares) are compared with exact quantum 
computation (circles). 



given momentum p (see Sec. As for the considered 
value of initial mean coordinate go = —0.2 the initial en- 
ergy is not very large (w 0.75Vb, where Vb is the height 
of the barrier) , the wave packet mainly stays in the well 
(for the time considered t = 400 only « 20% of the wave 
packet is transmitted, see Fig. ^ and the distribution in 
coordinate space is more compact. 

In Fig. [7|we present the dependence of tunneling time 
on initial mean position of the wave packet. Tunneling 
time is determined as the difference of presence times 
(EOI for points Xa = 0.5 • 0.6709 and Xb = 2.0 • 0.6709 (at 
X = 0.6709 potential has the maximum). Tunneling is 
usually stronger for the higher energy. Therefore, the in- 
crease of I go I (and corresponding increase of initial mean 
energy) leads to the growth of the average speed both of 
the transmitted part and of the wave packet as a whole. 
The transmitted part passes the region of the barrier (the 
space between the points Xa and Xb) faster, and so one 
expects that the increase of |go| causes the decrease of 
tunneling time. Indeed, the value of tunneling time drops 
with increase of |qo|- Results of QT simulation (squares 
in Fig. [TJ deviate from those of exact computation (cir- 
cles) within the range of errors. The deviation is maximal 
for large \qo\. Probably, this is because for large |<7o| the 
wave packet leaves the well almost entirely (see Fig. ^| , 
and the evolution of most trajectories corresponds to the 
unbounded accelerated motion. In such situation trajec- 
tories scatter and approximation H12|) does not represent 
the quantum tomogram as exactly as for smaller |go|- 



IV. SIMULATION OF THE EXCITON 
TUNNELING 

The method described in Sec. can be used to sim- 
ulate the evolution of systems with more than one de- 
gree of freedom. In this section we demonstrate such 
a possibility, considering nonstationary tunneling of the 
composite particle, exciton, through the potential barrier 
in one-dimensional (ID) semiconductor structure (quan- 
tum wire). Exciton is a bound state of electron and hole 
in semiconductor, therefore we deal with two degrees of 
freedom in contrast to Sec. IIIII 

Possible experimental realization is as follows. Con- 
sider semiconductor quasi-one-dimensional nanostruc- 
ture where the motion is allowed only in one direction 
(quantum wire). Transverse motion is restricted due to 
strong confining barriers. Potential barrier in the direc- 
tion of allowed motion can be located at some point of the 
quantum wire either using the semiconductor heterojunc- 
tion or by the gate. Using femtosecond laser pulses, we 
can form an excitonic wave packet, either by the quasires- 
onance pumping, or exciting an electron from valence 
band with formation of a hole and subsequent binding 
of two particles into exciton. Then the excitonic wave 
packet can move to the barrier and with the help of some 
detectors one can investigate scattering of the exciton. 

Keeping this in mind let us construct the model for 
the simulation. We use the constants corresponding to 
GaAs for reference (dielectric constant e — 12.5, effective 
masses of electron and hole are me = O.QlrrSe^ and = 
0.15mi°\ respectively, where mi"' is the electron mass 
in vacuum). 3D exciton in bulk GaAs is characterized 
by effective Bohr radius a* ~ lOnm and binding energy 
Ki AmeV. We use unit of length a* , unit of mass me, 
and h — 1. Corresponding units of energy and time are 
Eq = / {rrieO?) w lOmeV and to = mea/h « 100/s. 

The energy spectrum and wave functions of relative 
electron and hole motion in 3D exciton are analogous to 
those of hydrogen atom. But this is not the case for ID 
exciton. First, electron-hole effective interaction poten- 
tial in quasi-lD structure is not Coulomb. Indeed, if the 
exciton size in the direction of allowed motion is much 
greater than the width of the quantum wire (in the trans- 
verse direction), then the adiabatic approximation is ap- 
plicable and 3D interaction potential must be averaged 
over the transverse degrees of freedom. Resulting ID ef- 
fective potential substantially differs from Coulomb (see 
[sSi] for the discussion of similar model). Second, corre- 
sponding energy spectrum and wave functions of electron 
and hole relative motion also change in comparison with 
the hydrogen-like states. We choose the wave function of 
the exciton ground state in Gaussian form. 

Excitonic wave packet can be represented as a Gaus- 
sian wave packet in center-of-mass coordinates: 

g-rV(2cr) ^-(R~xof/{2S)+iRpo 

n-.,x,,t^o)^^-^ — ^5)v^ — ' (21) 
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where R = {rrieXe + mhXh)/{me + mji), r = \xe - Xh\, 
Xe and Xh are electron and hole coordinates, XqjPo and 
S are parameters; for them we used the following values: 
xo = — 10, po = 3, S' = 2 and a — 1. 

External potential is assumed to be zero everywhere 
except the region of barrier; we use the barrier of thick- 
ness equal to 5nm, or 0.5 in accepted units. For sim- 
plicity we set the barriers for electron and hole to be the 
same and use both external and interaction potentials 
in quadratic form, cut at some distance. Then external 
potential is given by 



Vext{x) 




if 



12^1 < 



if bl > J% 



(22) 



C is the height of the barrier, its width is ^JC jD — 0.5. 

Interaction potential Vint is also assumed to be 
quadratic: 




if 7- > , 



(23) 



where r = \X(^ — Xh\- Potential H23|) can describe, e.i?., 
e-h interaction in spatially indirect exciton, for example 
in coupled quantum wires with large interwire separation 
|36j | . Initial wave function of relative motion, chosen to be 
Gaussian with unity dispersion, is negligible within one 
percent accuracy at r = 3. Thus we ch oose the radius of 
electron- hole interaction to be yj A/B = 3. 

We assume that we deal with a quasi- ID exciton with 
binding energy Eq — 1/8. In fact, for an exciton in 
quantum wire, the wave function, binding energy, etc., 
are essentially influenced by the properties of quantum 
wire. We also neglect the possibility of electron and hole 
recombination at the time scales studied. 

Here we consider just an example, therefore use rela- 
tively simple model. Still this model contains the main 
features of exciton tunneling, such as the possibility of 
ionization, barrier and interaction with realistic strength 
and size, the fact that composite particle is the bounded 
state of two particles. 

For stationary state the binding energy is —Eq — 
I '^int{'^)Hint{r)^int{r)dr, where 5'i„t(r) is the wave 
function of relative motion. Then, from Ea. (|23|l and con- 
dition y/A/B = 3 we have A « ISEc/U. 

We do not use the variables of relative and center-of- 
mass motion, for QT simulation it is easier to deal with 
the initial conditions and evolution equation in coordi- 
nates of electron and hole x^ and x^. Potentials (|22|l and 
(123(1 are quadratic, that makes tomographic considera- 
tion of the problem easier (see Sec. IIII|) . discontinuity is 
neglected. On the other hand, in coordinates Xe and Xh 
the evolution equations depend on trajectory distribution 
(see Secnj, therefore the problem considered allows to 
employ all techniques, developed for one degree of free- 
dom, in this case of two degrees of freedom. 
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FIG. 8: Probability density distributions in coordinate space 
for electron {pe{x)) and hole {phix)) at times i = and t = 10. 
QT simulation (solid lines) is compared with exact numerical 
solution (dashed lines). All values are in units h = ml = 
Ec ~ 1, m* is electron effective mass and Ec is binding 
energy of the exciton. The height of the barrier C = 1, width 
^/C/D = 0.5. 
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FIG. 9: Probability density distributions in coordinate space 
for electron {pe{x)) and hole [phix)) at times t = 2Q and 
t = 30. QT simulation (solid lines) is compared with exact 
numerical solution (dashed lines) . The same units and barrier 
parameters as in Fig. Elare used. 



Results of exciton tunneling simulation are presented 
in Figs. 1811111 All parameters are fixed (see above), ex- 
cept the barrier height C in Fig.^l In Figs. |H1 and 121 we 
depict the evolution of probability density for electron 
and hole, barrier height C — \. Solid lines in Figs.|Hl- 
|51 and circles in Fig. 1101 correspond to the simulation in 
quantum tomography approach, while dashed lines and 
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FIG. 10: Probability of exciton ionization Pion due to elec- 
tron and hole scattering on the barrier in opposite directions 
versus barrier height C. Circles and squares represent QT 
simulation and exact solution, respectively. Considered is the 
barrier of thickness 0.5. The units are the same as in Fig. |H] 



squares represent the exact numerical computation. Un- 
like in Figs.[21-|ni here we show the results of several com- 
bined QT simulation runs, therefore corresponding lines 
are relatively smooth. Coincidence of QT and exact com- 
putation results, initially very good, becomes poorer with 
time (Figs-IHlandini), but QT simulation reproduces main 
properties of exciton tunneling: wave packets broadening 
with time and due to interaction with the barrier, shrink- 
ing near the barrier, dividing into two parts (reflected 
and transmitted). Note that QT results are quite close 
to exact ones even for long times {t = 30), despite the 
fact that the motion is unbounded (see Sec^J- Integral 
values (in Fig. I1U|I obtained in QT approach also agree 
with the exact results. Larger discrepancies correspond 
to higher barriers, probably due to larger inaccuracies, in- 
troduced by neglecting the potential discontinuity in the 
case of stronger interaction with the external potential. 

The electron and hole wave packets begin the motion 
from the point x = — 10 and, shrinking near the bar- 
rier, are partially reflected and transmitted. For the case 
presented in Figs. |H1- El about the half of wave pack- 
ets is transmitted. Interesting is the question about the 
ionization probability of exciton, induced by interaction 
with the barrier. If electron and hole are scattered in 
different directions on the barrier, the distance between 
them can become quite large, but, in principle, there is a 
possibility that exciton is not ionized after such scatter- 
ing, because one of the particles can be 'pulled' beyond 
the barrier, to the other particle, due to electron-hole 
attraction. On the other han d, the electron-hole interac- 
tion is cut at the distance y/ A/B in our model. After 
the interaction with the barrier the wave packet divides 
into reflected and transmitted parts moving in opposite 



directions. For the time large enough, these two parts 
are well separated, the separation between them grows 
and the leakage through the barrier in both directions 
is negligible. Denote the probability of ionization due 
to electron and hole scattering in different directions as 
Pion- Then, the probability to find electron and hole in 
different directions in respect to the barrier, with e-h dis- 
tance being larger than A/ B, approaches Pion m the 
limit t — > oo. 

The probability of ionization due to electron and hole 
scattering in different directions on the barrier Pjon is 
presented in Fig. 1101 depending on the barrier height C. 
For very high and very low barriers Pion must approach 
zero, because in former case both particles are reflected 
and in the latter they both are transmitted. This trend 
is seen in Fig. and Pion depending on C is max- 
imal (other parameters are fixed, see above) at C « 1. 
Note that these features are obvious for curves represent- 
ing both QT simulation (circles) and exact computation 
(squares), and in general two curves are quite close to 
each other. 



V. CONCLUSION 

We have developed the new method of numerical simu- 
lation of quantum nonstationary processes and applied it 
to the problem of tunneling of the wave packet through 
the potential barrier. The method is based on tomo- 
graphic representation of quantum mechanics. The quan- 
tum tomogram is used in a sense as the distribution 
function for the ensemble of trajectories in space X, /i, i/, 
where X = + vp is the coordinate measured in ro- 
tated and scaled reference frame, q,p are coordinate and 
momentum of the system, respectively. The trajectories 
are governed by the equations, resembling the Hamilton 
equations of motion, therefore, some analogue of molecu- 
lar dynamics can be used. The Gaussian approximation 
allows to avoid the direct calculation of quantum tomo- 
gram. Instead of quantum tomogram, the parameters of 
the approximation are used in the equations of motion. 
Those parameters can be obtained if one calculates the 
local moments of the ensemble of the trajectories. 

The problem of nonstationary tunneling of the wave 
packet was considered. Our method gave the results in 
agreement with those obtained by the method of " Wigner 
trajectories" and by exact quantum computation. 

Of course, we made only the first step in the develop- 
ment of this simulation method, having considered one- 
dimensional problem. But it is obvious that the general- 
ization for the case of several dimensions is straightfor- 
ward. In this case the method possesses the additional 
advantage of more rapid convergence due to the fact the 
quantum tomogram is nonnegative. This may help to 
overcome the problem of sign for fermionic systems. In 
the next work we intend to consider the many-body prob- 
lem for fermionic and bosonic systems by means of the 
quantum tomography method. The similar approach (in 
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the framework of "Wigner trajectories") has been already 
appUed for the investigation of the tunnehng of two iden- 
tical particles ■ 
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